#include <TROOT.h>
void p2FIT( char *momentum_order, int beta100, int l, int t, int nclust, char *dir, char *epm, char *Oa2print, char *rgiprint, int iprop, float p2min_m1, float p2max_m1, float p2min_m2, float p2max_m2, float cut)
{
        L=l;
        T=t;
        Nclust=nclust;

     float **p2_interval=(float**)malloc(sizeof(float*)*(2)); // M1 & M2
     for (int im=0; im<2; im++) p2_interval[im]=(float*)malloc(sizeof(float)*(2));  // begin & end

	p2_interval[0][0]=p2min_m1;
	p2_interval[0][1]=p2max_m1;
	p2_interval[1][0]=p2min_m2;
	p2_interval[1][1]=p2max_m2;

	for (int im=0; im<2; im++) printf("M%d %8.5f : %8.5f \n", im+1,p2_interval[im][0], p2_interval[im][1]); 
	printf("cut=%8.5f\n",cut);
   /////////////////////////
  //// Momenta set
  //////////////////////////////
        int set_p[4];
        char setpprint[16];

    
    whichProp(iprop,setpprint);

    char momentafile[1024]; sprintf(momentafile,"../INPUT/mom_%d.dat",beta100);
    read_momenta_data(momentafile);
    double *nn=(double*)malloc(sizeof(double)*(Nmom));
    compute_nn(momentum_order, nn);
    if (iprop<16)
    {
    whichSetp(iprop,set_p);
    Momentum_Magro(momentum_order, set_p);
    }
    else
    {
    Momentum(momentum_order);
    }
    // equivalent momenta
    Nmom_eq=Number_of_Equivalent_Momenta(SinP2);
    Average_over_Equivalent_Momenta(SinP2,SinP2,SinP2_eq);
    Average_over_Equivalent_Momenta(SinP2,SinP4,SinP4_eq);
    double *SinP4overSinP2Square_eq=malloc(sizeof(double)*(Nmom_eq));
     double *aux_mom=malloc(sizeof(double)*Nmom);
    for (int imom=0; imom<Nmom; imom++) aux_mom[imom]=SinP4[imom]/pow(SinP2[imom],2);
    Average_over_Equivalent_Momenta(SinP2,aux_mom,SinP4overSinP2Square_eq);

	//////////////////////////////////
	// IO files
	///////////////////////////////

	char file[1024];
	char gammaprint[8];
	char dir_qav[1024];

     	sprintf(dir_qav,"%s/%s/QAV/",dir,setpprint);

	//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//


	//////////////////////////////////
	// Read Zs
	/////////////////////////////
	printf("====================\n");
	printf(" READ Z             \n");

        double **Zq=malloc(sizeof(double*)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++) Zq[iclust]=malloc(sizeof(double)*(Nmom_eq));

        double ***ZO=malloc(sizeof(double**)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++) ZO[iclust]=malloc(sizeof(double*)*(Nmom_eq));
        for (int iclust=0; iclust<Nclust+1; iclust++)  for (int imom=0; imom<Nmom_eq; imom++) ZO[iclust][imom]=malloc(sizeof(double)*(5));
	double *aux_clust=malloc(sizeof(double)*(Nclust+1));


	// Zq

			for (int imom=0; imom<Nmom_eq; imom++)
			{
				sprintf(file,"%s/p%03d/Zq_%s%s.jk",dir_qav,imom,Oa2print,rgiprint);
				FILE *fin=fopen(file,"r");
				if(fin == (FILE *)NULL) {
					fprintf(stderr,"Unable to open file %s for reading\n",file);
				}

				for (int iclust=0; iclust<Nclust+1; iclust++) fscanf(fin,"%lf",&aux_clust[iclust]);
				fclose(fin);

				for (int iclust=0; iclust<Nclust+1; iclust++) Zq[iclust][imom]=aux_clust[iclust];
			}

	// ZO 
			for (int igamma=0; igamma<5; igamma++)
			{
				whichGammaProm(igamma,gammaprint);
				{
					for (int imom=0; imom<Nmom_eq; imom++)
					{

						sprintf(file,"%s/p%03d/Z%s_%s%s.jk",dir_qav,imom,gammaprint,Oa2print,rgiprint);
						FILE *fin=fopen(file,"r");
						if(fin == (FILE *)NULL) {
							fprintf(stderr,"Unable to open file %s for reading\n",file);
						}

						for (int iclust=0; iclust<Nclust+1; iclust++) fscanf(fin,"%lf",&aux_clust[iclust]);
						fclose(fin);

						for (int iclust=0; iclust<Nclust+1; iclust++) ZO[iclust][imom][igamma]=aux_clust[iclust];

					}

				}
			}


	printf("============================\n");
	printf(" p^2 fit  \n");
	// auxiliar for minuit routines
	// in practice I'll need Nmom_for_fit < Nmom 
	double *ZAUX=malloc(sizeof(double)*(Nmom_eq));
	double **ZAUX_CLUST=malloc(sizeof(double*)*(Nmom_eq));
	for (int imom=0; imom<Nmom_eq; imom++) ZAUX_CLUST[imom]=malloc(sizeof(double)*(Nclust+1));

	double *P2AUX=malloc(sizeof(double)*(Nmom_eq));

	double ***Z=malloc(sizeof(double**)*(6)); // One extra for Zq
	for (int igamma=0; igamma<6 ; igamma++)  Z[igamma]=malloc(sizeof(double*)*(2)); // M1 & M2
	for (int igamma=0; igamma<6 ; igamma++)  for (int im=0; im<2; im++) Z[igamma][im]=malloc(sizeof(double)*(Nclust+1));

	double ***slope=malloc(sizeof(double**)*(6)); // One extra for slope Zq
	for (int igamma=0; igamma<6 ; igamma++)  slope[igamma]=malloc(sizeof(double*)*(2)); // M1 & M2
	for (int igamma=0; igamma<6 ; igamma++)  for (int im=0; im<2; im++) slope[igamma][im]=malloc(sizeof(double)*(Nclust+1));

	double ***chi2=malloc(sizeof(double**)*(6)); // One extra for chi2 Zq
	for (int igamma=0; igamma<6 ; igamma++)  chi2[igamma]=malloc(sizeof(double*)*(2)); // M1 & M2
	for (int igamma=0; igamma<6 ; igamma++)  for (int im=0; im<2; im++) chi2[igamma][im]=malloc(sizeof(double)*(Nclust+1));


	int Nmom_for_fit;	


	///////////////////////////////////
	////////////////////////////
	// p^2 fit for Zq
	////////////////////////////////
	////////////////////////////////////

	printf("Zq\n");
		for (int im=0; im<2; im++)
		{
			printf("M%d\n",im+1);
			for (int iclust=0; iclust<Nclust+1; iclust++)
			{
				// Select momenta 
				int imom_for_fit=0;
				for (int imom=0; imom<Nmom_eq; imom++) 
				{
					if (SinP2_eq[imom] > p2_interval[im][0]  && SinP2_eq[imom]< p2_interval[im][1] && SinP4_eq[imom]/pow(SinP2_eq[imom],2)<cut )
//					if (SinP2_eq[imom] > p2_interval[im][0]  && SinP2_eq[imom]< p2_interval[im][1] && SinP4overSinP2Square_eq[imom]<cut)
					{
				//		if (iclust==0) printf("SinP2=%8.5f SinP4 %8.5f SinP4_eq/(SinP2_eq)^2=%8.5f vs [SinP4/(SinP2)^2]_eq=%8.5f\n",SinP2_eq[imom],SinP4_eq[imom],  SinP4_eq[imom]/pow(SinP2_eq[imom],2), SinP4overSinP2Square_eq[imom]);
						ZAUX[imom_for_fit]=Zq[iclust][imom];
						for (int iclust_loop=0; iclust_loop<Nclust+1; iclust_loop++) ZAUX_CLUST[imom_for_fit][iclust_loop]=Zq[iclust_loop][imom];
						P2AUX[imom_for_fit]=SinP2_eq[imom];
						imom_for_fit++;
					}
				}
				Nmom_for_fit=imom_for_fit;
/*                                        if (iclust==Nclust) {
                                                printf("Nmom_for_fit=%d\n",Nmom_for_fit);
                                                for (int imom=0; imom<Nmom_for_fit; imom++) printf("%8.5f %8.5f\n" , P2AUX[imom], ZAUX[imom]);
                                         }*/

				// M1
				if (im==0) lin_fit(ZAUX,ZAUX_CLUST,P2AUX,Nmom_for_fit,Nclust,Z[5][im][iclust],slope[5][im][iclust],chi2[5][im][iclust]); // Zq is saved in the last element (igamma=5)
				// M2
				//if (im==1) {wprom_fit(ZAUX,ZAUX_CLUST,P2AUX,Nmom_for_fit,Nclust,Z[5][im][iclust],chi2[5][im][iclust]); slope[5][im][iclust]=0.;}
				if (im==1) {wprom(ZAUX,ZAUX_CLUST,P2AUX,Nmom_for_fit,Nclust,Z[5][im][iclust]); slope[5][im][iclust]=0.;}
			}
                        sprintf(file,"%s/p2FIT/M%d/Zq_%s%s.fit",dir_qav,im+1,Oa2print,rgiprint);
                        FILE *fout=fopen(file,"w");
			for (int imom=0; imom<Nmom_for_fit; imom++)
                        {
				for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=ZAUX_CLUST[imom][iclust];
				fprintf(fout,"%8.5f %8.5f %8.5f\n",P2AUX[imom],ZAUX[imom],error_jack_real(aux_clust,Nclust));
			}
			for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=Z[5][im][iclust];
			fprintf(fout,"%8.5f %8.5f\n",Z[5][im][Nclust],error_jack_real(aux_clust,Nclust));
			for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=slope[5][im][iclust];
                        fprintf(fout,"%8.5f %8.5f\n",slope[5][im][Nclust],error_jack_real(aux_clust,Nclust));
			fprintf(fout,"%8.5f\n",chi2[5][im][Nclust]);
			fclose(fout);


		}

	/////////////////////////////////
	////////////////////////////
	// p^2 fit for ZO
	////////////////////////////////
	//////////////////////////////////


		for (int igamma=0; igamma<5; igamma++) {
		whichGammaProm(igamma,gammaprint);
		printf("Z%s\n",gammaprint);
			for (int im=0; im<2; im++){
				for (int iclust=0; iclust<Nclust+1; iclust++)
				{
					// Select momenta 
					int imom_for_fit=0;
					for (int imom=0; imom<Nmom_eq; imom++)
					{
						if (SinP2_eq[imom] > p2_interval[im][0]  && SinP2_eq[imom]< p2_interval[im][1] && SinP4_eq[imom]/pow(SinP2_eq[imom],2)<cut)
						{
							ZAUX[imom_for_fit]=ZO[iclust][imom][igamma];
							for (int iclust_loop=0; iclust_loop<Nclust+1; iclust_loop++) ZAUX_CLUST[imom_for_fit][iclust_loop]=ZO[iclust_loop][imom][igamma];
							P2AUX[imom_for_fit]=SinP2_eq[imom];
							imom_for_fit++;
						}

					}
					Nmom_for_fit=imom_for_fit;
	                /*                if (iclust==Nclust) {
                	                        printf("Nmom_for_fit=%d\n",Nmom_for_fit);
                        	                for (int imom=0; imom<Nmom_for_fit; imom++) printf("%8.5f %8.5f\n" , P2AUX[imom], ZAUX[imom]);
        	                         }
			*/
					// M1
					if (im==0) lin_fit(ZAUX,ZAUX_CLUST,P2AUX,Nmom_for_fit,Nclust,Z[igamma][im][iclust],slope[igamma][im][iclust],chi2[igamma][im][iclust]); // Zq is saved in the last element (igamma=igamma)
					// M2
					//if (im==1) {wprom_fit(ZAUX,ZAUX_CLUST,P2AUX,Nmom_for_fit,Nclust,Z[igamma][im][iclust],chi2[igamma][im][iclust]); slope[igamma][im][iclust]=0.;}
					if (im==1) {wprom(ZAUX,ZAUX_CLUST,P2AUX,Nmom_for_fit,Nclust,Z[igamma][im][iclust]); slope[igamma][im][iclust]=0.;}
				}
                        sprintf(file,"%s/p2FIT/M%d/Z%s_%s%s.fit",dir_qav,im+1,gammaprint,Oa2print,rgiprint);
                        FILE *fout=fopen(file,"w");
                        for (int imom=0; imom<Nmom_for_fit; imom++)
                        {
                                for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=ZAUX_CLUST[imom][iclust];
                                fprintf(fout,"%8.5f %8.5f %8.5f\n",P2AUX[imom],ZAUX[imom],error_jack_real(aux_clust,Nclust));
                        }
                        for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=Z[igamma][im][iclust];
                        fprintf(fout,"%8.5f %8.5f\n",Z[igamma][im][Nclust],error_jack_real(aux_clust,Nclust));
                        for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=slope[igamma][im][iclust];
                        fprintf(fout,"%8.5f %8.5f\n",slope[igamma][im][Nclust],error_jack_real(aux_clust,Nclust));
                        fprintf(fout,"%8.5f\n",chi2[igamma][im][Nclust]);
                        fclose(fout);

			}
		}

	// write output

	double **aux_allmomenta=malloc(sizeof(double*)*(Nmom_eq));
	for (int imom=0; imom<Nmom_eq; imom++) aux_allmomenta[imom]=malloc(sizeof(double)*(Nclust+1));

	printf("============================\n");
	printf("WRITING OUTPUT in %s/p2FIT/\n",dir_qav);    
	// Zq
	for (int im=0; im<2; im++)
	{
			for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=Z[5][im][iclust];
			sprintf(file,"%s/p2FIT/M%d/Zq_%s%s",dir_qav,im+1,Oa2print,rgiprint);
			write_output_fitmomenta_real(file,aux_clust,p2_interval[im][0], p2_interval[im][1],cut);

			for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=slope[5][im][iclust];
			sprintf(file,"%s/p2FIT/M%d/slope_Zq_%s%s",dir_qav,im+1,Oa2print,rgiprint);
			write_output_fitmomenta_real(file,aux_clust,p2_interval[im][0], p2_interval[im][1],cut);

	}



	// ZO
	for (int im=0; im<2; im++)
	{
			for (int igamma=0; igamma<5; igamma++)
			{
				whichGammaProm(igamma,gammaprint);
				for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=Z[igamma][im][iclust];
				sprintf(file,"%s/p2FIT/M%d/Z%s_%s%s",dir_qav,im+1,gammaprint,Oa2print,rgiprint);
				write_output_fitmomenta_real(file,aux_clust,p2_interval[im][0], p2_interval[im][1],cut);

				whichGammaProm(igamma,gammaprint);
				for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=slope[igamma][im][iclust];
				sprintf(file,"%s/p2FIT/M%d/slope_Z%s_%s%s",dir_qav,im+1,gammaprint,Oa2print,rgiprint);
				write_output_fitmomenta_real(file,aux_clust,p2_interval[im][0], p2_interval[im][1],cut);


			}
	}


	delete Zq; delete ZO;
	printf(" p2FIT.c finished succesfully! \n");

}
